/* Copyright 2012 Tobias Marschall
 *
 * This file is part of CLEVER.
 *
 * CLEVER is free software: you can redistribute it and/or modify
 * it under the terms of the GNU General Public License as published by
 * the Free Software Foundation, either version 3 of the License, or
 * (at your option) any later version.
 *
 * CLEVER is distributed in the hope that it will be useful,
 * but WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 * GNU General Public License for more details.
 *
 * You should have received a copy of the GNU General Public License
 * along with CLEVER.  If not, see <http://www.gnu.org/licenses/>.
 */

#include <sstream>
#include <math.h>

#include "Read.h"
#include "SplitMapping.h"
#include "../BamHelper.h"

using namespace std;

SplitMapping::SplitMapping(Read& parent) : left_anchor_idx(-1), right_anchor_idx(-1), parent(parent), alignments(0) {
}

SplitMapping::SplitMapping(int left_anchor_idx, int right_anchor_idx, Read& parent, vector<SplitAligner::split_alignment_t>* alignments, int phred_score, int mismatch_phred_score) : Mapping(phred_score, mismatch_phred_score), left_anchor_idx(left_anchor_idx), right_anchor_idx(right_anchor_idx), parent(parent), alignments(alignments) {
	if ((alignments != 0) && (alignments->size() > 0)) {
		phred_score = alignments->at(0).phred_score;
		assert(phred_score >= 0);
		probability = pow(10.0, -((double)phred_score)/10.0);
	}
}

SplitMapping::~SplitMapping() {
	if (alignments != 0) delete alignments;
}

BamTools::BamAlignment SplitMapping::getBamAlignment(Read& parent, const string& name) const {
	assert(alignments != 0);
	assert(alignments->size() > 0);
	const SplitAligner::split_alignment_t& split_aln = alignments->at(0);
	assert(left_anchor_idx < (int)parent.leftAnchors().size());
	const Read::anchor_t& left_anchor = parent.leftAnchors()[left_anchor_idx];
	BamTools::BamAlignment a;
	a.Name = name;
	a.Position = left_anchor.position;
	a.RefID = left_anchor.ref_id;
	if (left_anchor.reverse) {
		a.SetIsReverseStrand(true);
		ShortDnaSequence rev = parent.getSequence().reverseComplement();
		a.QueryBases = rev.toString();
		a.Qualities = rev.qualityString();
	} else {
		a.SetIsReverseStrand(false);
		a.QueryBases = parent.getSequence().toString();
		a.Qualities = parent.getSequence().qualityString();
	}
	a.MapQuality = 255;
	a.CigarData = split_aln.cigar;
	a.SetIsPaired(true);
	assert(phred_score >= 0);
	uint32_t as_tag = phred_score;
	if (!a.AddTag("AS", "I", as_tag)) assert(false);
	uint32_t ym_tag = mismatch_phred_score;
	if (!a.AddTag("YM", "I", ym_tag)) assert(false);
	// include YA tag for alternative CIGAR strings
	if (alignments->size() > 1) {
		ostringstream oss;
		for (size_t i = 1; i < alignments->size(); ++i) {
			if (i>1) oss << ';';
			oss << alignments->at(i).phred_score << ',' << alignments->at(i).mismatch_phred_score << ',' << alignments->at(i).cigar;
		}
		if (!a.AddTag("YA", "Z", oss.str())) assert(false);
	}
	return a;
}

int SplitMapping::getRefId() const {
	assert(left_anchor_idx < (int)parent.leftAnchors().size());
	return parent.leftAnchors()[left_anchor_idx].ref_id;
}

int SplitMapping::getStartPosition() const {
	assert(left_anchor_idx < (int)parent.leftAnchors().size());
	return parent.leftAnchors()[left_anchor_idx].position;
}

int SplitMapping::getEndPosition() const {
	assert(right_anchor_idx < (int)parent.rightAnchors().size());
	return parent.rightAnchors()[right_anchor_idx].position;
}

bool SplitMapping::isReverse() const {
	assert(left_anchor_idx < (int)parent.leftAnchors().size());
	return parent.leftAnchors()[left_anchor_idx].reverse;
}

int SplitMapping::longestIndel() const {
	assert(alignments != 0);
	assert(alignments->size() > 0);
	int result = 0;
	const vector<BamTools::CigarOp>& cigar = alignments->at(0).cigar;
	vector<BamTools::CigarOp>::const_iterator it = cigar.begin();
	for (; it != cigar.end(); ++it) {
		if ((it->Type == 'I') || (it->Type == 'D')) {
			result = max(result, (int)it->Length);
		}
	}
	return result;
}

auto_ptr<vector<Variation> > SplitMapping::getPutativeVariations(const std::vector<NamedDnaSequence*>& reference_list, bool rightmost) const {
	assert(alignments != 0);
	assert(alignments->size() > 0);
	auto_ptr<vector<Variation> > result(new vector<Variation>());
	const SplitAligner::split_alignment_t& aln = alignments->at(0);
	if (aln.type == SplitAligner::SPLIT_DELETION) {
		result->push_back(Variation(reference_list[getRefId()]->getName(), aln.reference_breakpoint - aln.length, aln.reference_breakpoint, -1.0, Variation::DELETION));
	}
	if (aln.type == SplitAligner::SPLIT_INSERTION) {
		string read_seq = isReverse()?parent.getSequence().reverseComplement().toString():parent.getSequence().toString();
		result->push_back(Variation(reference_list[getRefId()]->getName(), aln.reference_breakpoint, aln.length, -1.0, read_seq.substr(aln.query_breakpoint - 1, aln.length), Variation::INSERTION));
	}
	assert(result->size() == 1);
// 	cerr << "Before canonification: " << result->at(0) << endl;
	result->at(0).canonify(reference_list[getRefId()], rightmost);
// 	cerr << "After canonification:  " << result->at(0) << endl;
	return result;
}

void SplitMapping::print(std::ostream& os) const {
	os << getRefId() << ", " << getStartPosition() << " - " << getEndPosition() << ", PHRED: " << getPhredScore() << ", " << (isReverse()?"REVCOMP":"FORWARD") << ", SPLIT";
}
